Curvas de supervivencia con R

Autor/a

Santiago García Sánchez

Introducción

Se comparan distintos paquetes para graficar curvas de supervivencia con R.

Los conjuntos de datos que se emplean para los ejemplos son los siguientes:

veteran (incluido en el paquete survival): contrasta los resultados de dos regímenes de tratamiento para el cáncer de pulmón.
bladder (paquete survival): contiene datos sobre la recurrencia del cáncer de vejiga.
pharmacoSmoking (paquete asaur): datos provenientes de un ensayo aleatorizado que compara la eficacia de la triple terapia frente al uso del parche para dejar de fumar.

Código

library(survival)
library(asaur)

veteran <- survival::veteran
bladder <- survival::bladder
pharmacoSmoking <- asaur::pharmacoSmoking

Con las funciones Surv(), survfit() y coxph() del paquete survival se ajustan distintas curvas de supervivencia:

Para el dataset veteran, se estiman funciones de supervivencia según tratamiento utilizando el método de Kaplan-Meier.
Para el dataset bladder, se obtiene la función de supervivencia con el estimador de Fleming-Harrington. Esta estimación será utilizada posteriormente para construir la función hazard estimada por el método de Nelson-Aalen.
Con el dataset pharmacoSmoking, se ajusta un modelo de hazards proporcionales de Cox considerando las variables grp (grupo asignado), ageGroup4 (edad) y employment (tipo de empleo). Luego, se obtiene la curva de supervivencia estimada por el modelo.

Código

km_fit <- survfit(Surv(veteran$time, veteran$status == 1) ~ veteran$trt,
                  type = "kaplan-meier", conf.int = 0.95, conf.type = "log-log",
                  error = "greenwood")

na_fit <- survfit(Surv(bladder$stop, bladder$event == 1) ~ 1,
                  type = "fh", conf.type = "log-log")

cox <- coxph(Surv(ttr, relapse) ~ grp + ageGroup4 + employment,
             data = pharmacoSmoking)
cox_fit <- survfit(cox)

R base y survival

Se realizan gráficos con R base y survival:

Código

plot(km_fit,
     main = "Función de supervivencia estimada por el método de Kaplan-Meier",
     xlab = "Tiempo (días)", ylab = "Probabilidad de supervivencia estimada",
     cex.main = 0.9, lty = 1, lwd = 2, cex = 1, cex.lab = 0.95, mark = 1,
     conf.int = TRUE, col = c("steelblue", "pink2"))

Código

plot(na_fit$time, main = "Función Hazard estimada por Nelson-Aalen",
     y = na_fit$n.event / na_fit$n.risk,
     type = "p", pch = 20, xlab = "Tiempo (días)", ylab = "Tasa")
lines(lowess(na_fit$time, na_fit$n.event / na_fit$n.risk), col = "violet")

Código

plot(cox_fit, main = "Modelo de regresión de Cox",
     cex.main = 0.9, lty = 1, lwd = 2, cex = 1, cex.lab = 0.95, mark = 1,
     col = "red3",
     ylab = "Probabilidad de supervivencia estimada", xlab = "Tiempo (días)")

Paquetes ggplot2 y ggfortify

Se utiliza la función autoplot.survfit() del paquete ggfortify y varias funciones del paquete ggplot2:

Código

library(ggplot2)
library(ggfortify)

autoplot(km_fit, conf.int.alpha = 0.2, censor.shape = "x") +
  labs(x = "Tiempo (días)",
       y = "Probabilidad de supervivencia estimada",
       title = "Función de supervivencia estimada por el método de Kaplan-Meier",
       color = "Tratamiento", fill = "Tratamiento") +
  theme(title = element_text(size = 10.5),
        plot.title = element_text(size = 11.5, face = "bold"))

Código

autoplot(cox_fit, conf.int.alpha = 0.2, conf.int.linetype = "dashed",
         conf.int.colour = "red", surv.colour = "red3", surv.size = 1.2,
         censor = FALSE)  +
  labs(x = "Tiempo (días)",
       y = "Probabilidad de supervivencia estimada",
       title = "Modelo de regresión de Cox")  +
  theme(title = element_text(size = 10.5),
        plot.title = element_text(size = 11.5, face = "bold"))

Paquete GGally

Se utiliza la función ggsurv() del paquete GGally, que genera curvas de supervivencia automáticamente:

Código

library(GGally)

ggsurv(km_fit)  +
  labs(x = "Tiempo (días)", y = "Probabilidad de supervivencia estimada",
       title = "Función de supervivencia estimada por el método de Kaplan-Meier",
       color = "Tratamiento", linetype = "Tratamiento")

Código

ggsurv(cox_fit, surv.col = "red3")  +
  labs(x = "Tiempo (días)",
       y = "Probabilidad de supervivencia estimada",
       title = "Modelo de regresión de Cox")

Paquete survminer

El paquete survminer ofrece la función ggsurvplot() para dibujar curvas de supervivencia:

Código

library(survminer)

ggsurvplot(km_fit, data = veteran, risk.table = TRUE, conf.int = TRUE,
           break.time.by = 50, xlim = c(0, 600), surv.median.line = "hv",
           ggtheme = theme_minimal())

Código

ggsurvplot(cox_fit, data = pharmacoSmoking, palette = "red3", risk.table = TRUE,
           conf.int = TRUE, surv.median.line = "hv", break.time.by = 50,
           ggtheme = theme_minimal())

Paquete ggsurvfit

Se construyen curvas de supervivencia con tablas de riesgo recurriendo a funciones del paquete ggsurvfit (principalmente ggsurvfit()).

Código

library(ggsurvfit)

survfit2(Surv(time, status) ~ trt, data = veteran) %>%
  ggsurvfit() +
  add_confidence_interval() +
  add_risktable(risktable_group = "risktable_stats") +
  scale_ggsurvfit() +
  add_risktable_strata_symbol(symbol = "\U25CF", size = 10)

Código

coxph(Surv(ttr, relapse) ~ grp + ageGroup4 + employment, data = pharmacoSmoking) %>%
  survfit2() %>%
  ggsurvfit() +
  add_confidence_interval() +
  add_risktable() +
  scale_ggsurvfit()  +
  scale_y_continuous(limits = c(0, 1))

Paquete rms

El paquete rms incluye las funciones psm() y npsurv() para modelos de supervivencia paramétricos y no paramétricos respectivamente. Las curvas ajustadas pueden graficarse con survplot() (visualizaciones estáticas) o survplotp() (visualizaciones interactivas).

Código

library(rms)

km_fit_rms <- npsurv(Surv(time, status == 1) ~ trt, data = veteran)

survplot(km_fit_rms,
         col = c("steelblue", "red"),
         xlab = "Tiempo (días)",
         ylab = "Probabilidad de supervivencia estimada")

Código

survplotp(km_fit_rms,
          col = c("steelblue", "pink2"),
          xlab = "Tiempo (días)",
          ylab = "Probabilidad de supervivencia estimada")

Código

survplot(psm(Surv(bladder$stop, bladder$event == 1) ~ 1),
         what = "hazard", xlab = "Tiempo (días)")

Paquete jskm

La función jskm() del paquete del mismo nombre construye un gráfico Kaplan-Meier junto con una tabla de riesgo.

Código

library(jskm)

jskm(km_fit, table = TRUE, ci = TRUE,
     main = "Función de supervivencia estimada por el método de Kaplan-Meier",
     ystrataname = "Tratamiento", ystratalabs = c("1", "2"),
     xlabs = "Tiempo (días)", ylabs = "Probabilidad de supervivencia estimada",
     label.nrisk = "N.º en riesgo")

Paquete KMunicate

El paquete KMunicate crea gráficos Kaplan-Meier utilizando un estilo que sigue una serie de recomendaciones:

Código

library(KMunicate)

KMunicate(km_fit, seq(0, max(veteran$time), length.out = 10))

--- title: "Curvas de supervivencia con R" author: "Santiago García Sánchez" format: html: code-fold: show code-tools: true code-link: true df-print: paged toc: true toc-location: left theme: spacelab editor: source lang: es self-contained: true execute: warning: false error: false --- ## Introducción Se comparan distintos paquetes para graficar curvas de supervivencia con R. Los conjuntos de datos que se emplean para los ejemplos son los siguientes: - **[veteran](https://rdrr.io/cran/survival/man/veteran.html)** (incluido en el paquete survival): contrasta los resultados de dos regímenes de tratamiento para el cáncer de pulmón. - **[bladder](https://rdrr.io/cran/survival/man/bladder.html)** (paquete survival): contiene datos sobre la recurrencia del cáncer de vejiga. - **[pharmacoSmoking](https://search.r-project.org/CRAN/refmans/asaur/html/pharmacoSmoking.html)** (paquete [asaur](https://cran.r-project.org/web/packages/asaur/index.html)): datos provenientes de un ensayo aleatorizado que compara la eficacia de la triple terapia frente al uso del parche para dejar de fumar. ```{r} library(survival) library(asaur) veteran <- survival::veteran bladder <- survival::bladder pharmacoSmoking <- asaur::pharmacoSmoking ``` Con las funciones `Surv()`, `survfit()` y `coxph()` del paquete [**survival**](https://github.com/therneau/survival) se ajustan distintas curvas de supervivencia: - Para el dataset `veteran`, se estiman funciones de supervivencia según tratamiento utilizando el método de Kaplan-Meier. - Para el dataset `bladder`, se obtiene la función de supervivencia con el estimador de Fleming-Harrington. Esta estimación será utilizada posteriormente para construir la función hazard estimada por el método de Nelson-Aalen. - Con el dataset `pharmacoSmoking`, se ajusta un modelo de hazards proporcionales de Cox considerando las variables `grp` (grupo asignado), `ageGroup4` (edad) y `employment` (tipo de empleo). Luego, se obtiene la curva de supervivencia estimada por el modelo. ```{r} km_fit <- survfit(Surv(veteran$time, veteran$status == 1) ~ veteran$trt, type = "kaplan-meier", conf.int = 0.95, conf.type = "log-log", error = "greenwood") na_fit <- survfit(Surv(bladder$stop, bladder$event == 1) ~ 1, type = "fh", conf.type = "log-log") cox <- coxph(Surv(ttr, relapse) ~ grp + ageGroup4 + employment, data = pharmacoSmoking) cox_fit <- survfit(cox) ``` ## R base y survival Se realizan gráficos con R base y survival: ::: panel-tabset ## Kaplan-Meier ```{r} plot(km_fit, main = "Función de supervivencia estimada por el método de Kaplan-Meier", xlab = "Tiempo (días)", ylab = "Probabilidad de supervivencia estimada", cex.main = 0.9, lty = 1, lwd = 2, cex = 1, cex.lab = 0.95, mark = 1, conf.int = TRUE, col = c("steelblue", "pink2")) ``` ## Nelson-Aalen ```{r} plot(na_fit$time, main = "Función Hazard estimada por Nelson-Aalen", y = na_fit$n.event / na_fit$n.risk, type = "p", pch = 20, xlab = "Tiempo (días)", ylab = "Tasa") lines(lowess(na_fit$time, na_fit$n.event / na_fit$n.risk), col = "violet") ``` ## Cox ```{r} plot(cox_fit, main = "Modelo de regresión de Cox", cex.main = 0.9, lty = 1, lwd = 2, cex = 1, cex.lab = 0.95, mark = 1, col = "red3", ylab = "Probabilidad de supervivencia estimada", xlab = "Tiempo (días)") ``` ::: ## Paquetes ggplot2 y ggfortify Se utiliza la función `autoplot.survfit()` del paquete [**ggfortify**](https://github.com/sinhrks/ggfortify) y varias funciones del paquete [**ggplot2**](https://ggplot2.tidyverse.org/): ::: panel-tabset ## Kaplan-Meier ```{r} library(ggplot2) library(ggfortify) autoplot(km_fit, conf.int.alpha = 0.2, censor.shape = "x") + labs(x = "Tiempo (días)", y = "Probabilidad de supervivencia estimada", title = "Función de supervivencia estimada por el método de Kaplan-Meier", color = "Tratamiento", fill = "Tratamiento") + theme(title = element_text(size = 10.5), plot.title = element_text(size = 11.5, face = "bold")) ``` ## Cox ```{r} autoplot(cox_fit, conf.int.alpha = 0.2, conf.int.linetype = "dashed", conf.int.colour = "red", surv.colour = "red3", surv.size = 1.2, censor = FALSE) + labs(x = "Tiempo (días)", y = "Probabilidad de supervivencia estimada", title = "Modelo de regresión de Cox") + theme(title = element_text(size = 10.5), plot.title = element_text(size = 11.5, face = "bold")) ``` ::: ## Paquete GGally Se utiliza la función `ggsurv()` del paquete [**GGally**](https://ggobi.github.io/ggally/index.html), que genera curvas de supervivencia automáticamente: ::: panel-tabset ## Kaplan-Meier ```{r, warning=FALSE} library(GGally) ggsurv(km_fit) + labs(x = "Tiempo (días)", y = "Probabilidad de supervivencia estimada", title = "Función de supervivencia estimada por el método de Kaplan-Meier", color = "Tratamiento", linetype = "Tratamiento") ``` ## Cox ```{r} ggsurv(cox_fit, surv.col = "red3") + labs(x = "Tiempo (días)", y = "Probabilidad de supervivencia estimada", title = "Modelo de regresión de Cox") ``` ::: ## Paquete survminer El paquete [**survminer**](https://rpkgs.datanovia.com/survminer/index.html) ofrece la función `ggsurvplot()` para dibujar curvas de supervivencia: ::: panel-tabset ## Kaplan-Meier ```{r, warning=FALSE, message=FALSE} library(survminer) ggsurvplot(km_fit, data = veteran, risk.table = TRUE, conf.int = TRUE, break.time.by = 50, xlim = c(0, 600), surv.median.line = "hv", ggtheme = theme_minimal()) ``` ## Cox ```{r} ggsurvplot(cox_fit, data = pharmacoSmoking, palette = "red3", risk.table = TRUE, conf.int = TRUE, surv.median.line = "hv", break.time.by = 50, ggtheme = theme_minimal()) ``` ::: ## Paquete ggsurvfit Se construyen curvas de supervivencia con tablas de riesgo recurriendo a funciones del paquete [**ggsurvfit**](https://www.danieldsjoberg.com/ggsurvfit/) (principalmente `ggsurvfit()`). ::: panel-tabset ## Kaplan-Meier ```{r, warning=F, message=F} library(ggsurvfit) survfit2(Surv(time, status) ~ trt, data = veteran) %>% ggsurvfit() + add_confidence_interval() + add_risktable(risktable_group = "risktable_stats") + scale_ggsurvfit() + add_risktable_strata_symbol(symbol = "\U25CF", size = 10) ``` ## Cox ```{r, warning=F, message=F} coxph(Surv(ttr, relapse) ~ grp + ageGroup4 + employment, data = pharmacoSmoking) %>% survfit2() %>% ggsurvfit() + add_confidence_interval() + add_risktable() + scale_ggsurvfit() + scale_y_continuous(limits = c(0, 1)) ``` ::: ## Paquete rms El paquete **[rms](https://github.com/harrelfe/rms)** incluye las funciones `psm()` y `npsurv()` para modelos de supervivencia paramétricos y no paramétricos respectivamente. Las curvas ajustadas pueden graficarse con `survplot()` (visualizaciones estáticas) o `survplotp()` (visualizaciones interactivas). ::: panel-tabset ## Kaplan-Meier ```{r} library(rms) km_fit_rms <- npsurv(Surv(time, status == 1) ~ trt, data = veteran) survplot(km_fit_rms, col = c("steelblue", "red"), xlab = "Tiempo (días)", ylab = "Probabilidad de supervivencia estimada") ``` ## Kaplan-Meier (interactivo) ```{r} survplotp(km_fit_rms, col = c("steelblue", "pink2"), xlab = "Tiempo (días)", ylab = "Probabilidad de supervivencia estimada") ``` ## Función hazard ```{r} survplot(psm(Surv(bladder$stop, bladder$event == 1) ~ 1), what = "hazard", xlab = "Tiempo (días)") ``` ::: ## Paquete jskm La función `jskm()` del [paquete del mismo nombre](https://jinseob2kim.github.io/jskm/) construye un gráfico Kaplan-Meier junto con una tabla de riesgo. ```{r} library(jskm) jskm(km_fit, table = TRUE, ci = TRUE, main = "Función de supervivencia estimada por el método de Kaplan-Meier", ystrataname = "Tratamiento", ystratalabs = c("1", "2"), xlabs = "Tiempo (días)", ylabs = "Probabilidad de supervivencia estimada", label.nrisk = "N.º en riesgo") ``` ## Paquete KMunicate El paquete [**KMunicate**](https://ellessenne.github.io/KMunicate-package/) crea gráficos Kaplan-Meier utilizando un estilo que sigue una [serie de recomendaciones](https://bmjopen.bmj.com/content/9/9/e030215): ```{r} library(KMunicate) KMunicate(km_fit, seq(0, max(veteran$time), length.out = 10)) ```